# Replication script for: "Attribution and Migration: A Survey Experiment"

# Authors: Faure, M., Kantorowicz, J., & Weiss, A

# Load necessary packages

library(tidyverse)
library(ggthemes)
library(rstatix)
library(ggridges)
library(ggstatsplot)
library(ggpubr)
library(ggsci)
library(stargazer)
library(effsize)
library(modelsummary)
library(kableExtra)
library(vtable)

# Load the datasets

experiment_main_data <- readRDS("data/experiment_main_data.RDS")

experiment_manipulation_check <- readRDS("data/experiment_manipulation_check.RDS")

## Perform balance tests

names(experiment_main_data)

### Age

age <- lm(formula = age ~ experimental_groups_a, data = experiment_main_data)

### Gender

female <- lm(formula = female ~ experimental_groups_a, data = experiment_main_data)

### Political ideology

political_pref <- lm(formula = pref_political_1 ~ experimental_groups_a, data = experiment_main_data)

### Climate change a serious problem

climate_change <- lm(formula = pref_climate_1_1 ~ experimental_groups_a, data = experiment_main_data)

### Education

education <- lm(formula = higher_education ~ experimental_groups_a, data = experiment_main_data)

### Put all regressions together: Table 2

stargazer(age, female, education, political_pref, climate_change,
          type = "html",
          title = "Regression Table",
          style = "qje",
          out = "tables_figures/balancing_test.doc",
          dep.var.labels = c("Age", "Female", "Higher education", "Political ideology", "Climate change attitudes"),
          covariate.labels = c("Blame NL", "Blame EU", "Blame US"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

## Descriptive statistics

summary_stats <- experiment_main_data %>% select(age, female, higher_education, pref_political_1,
                                               pref_climate_1_1, outcome_blame)

summary_stats <- summary_stats %>% rename(Age = age,
                                          Female = female,
                                          `Higher education` = higher_education,
                                          `Political preferences` = pref_political_1,
                                          `Climate change attitudes` = pref_climate_1_1,
                                          `Favorability` = outcome_blame)                                       

### Put all regressions together: Table 1

st(data = summary_stats,
   out = "csv",
   file = "tables_figures/summary_stats_factors.csv",
   add.median = TRUE)

## Main figure

experiment_main_data$experimental_groups_a <- factor(experiment_main_data$experimental_groups_a,
                                                   levels = c("Control",
                                                              "Blame NL",
                                                              "Blame EU",
                                                              "Blame US"))

experiment_manipulation_check$experimental_groups_a <- factor(experiment_manipulation_check$experimental_groups_a,
                                                              levels = c("Control",
                                                                         "Blame NL",
                                                                         "Blame EU",
                                                                         "Blame US"))

p <- ggboxplot(data = experiment_main_data, 
               x = "experimental_groups_a", 
               y = "outcome_blame",
               color = "experimental_groups_a", 
               palette = get_palette("aaas", 4)) +
  geom_jitter(alpha=0.1, width = 0.3, height = 0.3) +
  theme(legend.position = "none") + 
  stat_compare_means(method = "t.test",
                     aes(label = ..p..),
                     comparisons = list(c("Blame NL", "Control"), 
                                        c("Blame NL", "Blame EU"),
                                        c("Blame NL", "Blame US"),
                                        c("Blame EU", "Blame US")),
                     label = "p.signif", symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("****", "***", "**", "*", "ns"))) +
  labs(x = "", 
       y = "Favorability") +
  scale_y_continuous(breaks=c(1,3,5,7)) +
  stat_summary(geom="errorbar", fun.data=mean_cl_normal, width=0.1, conf.int=0.95, color="black") +
  stat_summary(geom="point", fun.y=mean, color="black", size=2) +
  ggtitle("A. Full sample") +
  theme(plot.title = element_text(face="bold", hjust = 0.5))

g <- ggboxplot(data = experiment_manipulation_check, 
               x = "experimental_groups_a", 
               y = "outcome_blame",
               color = "experimental_groups_a", 
               palette = get_palette("aaas", 4)) +
  geom_jitter(alpha=0.1, width = 0.3, height = 0.3) +
  theme(legend.position = "none") + 
  stat_compare_means(method = "t.test",
                     aes(label = ..p..),
                     comparisons = list(c("Blame NL", "Control"), 
                                        c("Blame NL", "Blame EU"),
                                        c("Blame NL", "Blame US"),
                                        c("Blame EU", "Blame US")),
                     label = "p.format", symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
                                                            symbols = c("****", "***", "**", "*", "ns"))) +
  labs(x = "", 
       y = "Favorability") +
  scale_y_continuous(breaks=c(1,3,5,7)) +
  stat_summary(geom="errorbar", fun.data=mean_cl_normal, width=0.1, conf.int=0.95, color="black") +
  stat_summary(geom="point", fun.y=mean, color="black", size=2) +
  ggtitle("B. Excluding manipulation \n check failures") +
  theme(plot.title = element_text(face="bold", hjust = 0.5))

figure <- ggarrange(p, g, ncol = 2, widths = c(1, 1))

### Figure 1

ggsave(figure, file = "tables_figures/main_attribution.jpg", 
       scale = 0.8, 
       width = 12, 
       height = 8)

## First regression

experiment_main_data$experimental_groups_a <- factor(experiment_main_data$experimental_groups_a,
                                                   levels = c("Blame NL", 
                                                              "Control", 
                                                              "Blame EU", 
                                                              "Blame US"))

experiment_main_data$pref_political_1 <- as.numeric(experiment_main_data$pref_political_1)

experiment_manipulation_check$experimental_groups_a <- factor(experiment_manipulation_check$experimental_groups_a,
                                                        levels = c("Blame NL", 
                                                                   "Control", 
                                                                   "Blame EU", 
                                                                   "Blame US"))

experiment_manipulation_check$pref_political_1 <- as.numeric(experiment_manipulation_check$pref_political_1)

model1 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_main_data)

model2 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_main_data,
             weights = weight)

model3 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_main_data)

model4 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_main_data,
             weights = weight)

model5 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_manipulation_check)

model6 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_manipulation_check)

### Put all regressions together: Table 3

stargazer(model1, model2, model3, model4, model5, model6,
          type = "html",
          title = "Regression Table",
          style = "qje",
          out = "tables_figures/regression_table_attribution.doc",
          dep.var.labels = c("Favorability towards migrants"),
          covariate.labels = c("Control", "Blame EU", "Blame US", "Age", "Female","Higher education", "Political preferences", "Climate change attitudes"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

## Second regression

experiment_main_data$experimental_groups_a <- factor(experiment_main_data$experimental_groups_a,
                                                   levels = c("Blame EU",
                                                              "Control",
                                                              "Blame NL",
                                                              "Blame US"))


experiment_manipulation_check$experimental_groups_a <- factor(experiment_manipulation_check$experimental_groups_a,
                                                        levels = c("Blame EU",
                                                                   "Control",
                                                                   "Blame NL",
                                                                   "Blame US"))

model1 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_main_data)

model2 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_main_data,
             weights = weight)

model3 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_main_data)

model4 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_main_data,
             weights = weight)

model5 <- lm(formula = outcome_blame ~ experimental_groups_a, 
             data = experiment_manipulation_check)

model6 <- lm(formula = outcome_blame ~ experimental_groups_a + age + gender_cat + higher_education + pref_political_1 + pref_climate_1_1, 
             data = experiment_manipulation_check)

### Put all regressions together: Table 4

stargazer(model1, model2, model3, model4, model5, model6,
          type = "html",
          title = "Regression Table",
          style = "qje",
          out = "tables_figures/regression_table_attribution_2.doc",
          dep.var.labels = c("Favorability towards migrants"),
          covariate.labels = c("Control", "Blame NL", "Blame US", "Age", "Female","Higher education", "Political preferences", "Climate change attitudes"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001))

## Mann-Whitney tests of differences and effect sizes

### Control vs NL: full sample

dataset_experiment_contr_nl <- experiment_main_data %>% 
  filter(experimental_groups_a=="Control" | 
           experimental_groups_a=="Blame NL")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_contr_nl)

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_contr_nl,
            weights = dataset_experiment_contr_nl$weight)

dataset_experiment_contr_nl$experimental_groups_a <- as.numeric(dataset_experiment_contr_nl$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_contr_nl)

### Control vs NL: manipulation

dataset_experiment_contr_nl <- experiment_manipulation_check %>% 
  filter(experimental_groups_a=="Control" | 
           experimental_groups_a=="Blame NL")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_contr_nl) 

dataset_experiment_contr_nl$experimental_groups_a <- as.numeric(dataset_experiment_contr_nl$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_contr_nl)

### US vs NL: full sample

dataset_experiment_us_nl <- experiment_main_data %>% 
  filter(experimental_groups_a=="Blame US" | 
           experimental_groups_a=="Blame NL")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_us_nl)

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_us_nl,
            weights = dataset_experiment_us_nl$weight)

dataset_experiment_us_nl$experimental_groups_a <- as.numeric(dataset_experiment_us_nl$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_nl)

### US vs NL: manipulation

dataset_experiment_us_nl <- experiment_manipulation_check %>% 
  filter(experimental_groups_a=="Blame US" | 
           experimental_groups_a=="Blame NL")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_us_nl) 

dataset_experiment_us_nl$experimental_groups_a <- as.numeric(dataset_experiment_us_nl$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_nl)

cohens_d(outcome_blame ~ experimental_groups_a,
         data = dataset_experiment_us_nl)

### US vs EU: full sample

dataset_experiment_us_eu <- experiment_main_data %>% 
  filter(experimental_groups_a=="Blame US" | 
           experimental_groups_a=="Blame EU")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_us_eu)

dataset_experiment_us_eu$experimental_groups_a <- as.numeric(dataset_experiment_us_eu$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_eu)

### Control vs EU: full sample

dataset_experiment_control_eu <- experiment_main_data %>% 
  filter(experimental_groups_a=="Control" | 
           experimental_groups_a=="Blame EU")

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_control_eu)

dataset_experiment_us_eu$experimental_groups_a <- as.numeric(dataset_experiment_us_eu$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_eu)

### US vs EU: manipulation

dataset_experiment_us_eu <- experiment_manipulation_check %>% 
  filter(experimental_groups_a=="Blame US" | 
           experimental_groups_a=="Blame EU")

dataset_experiment_us_eu$experimental_groups_a <- as.numeric(dataset_experiment_us_eu$experimental_groups_a)
dataset_experiment_us_eu$outcome_blame <- as.numeric(dataset_experiment_us_eu$outcome_blame)

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_us_eu) 

dataset_experiment_us_eu$experimental_groups_a <- as.numeric(dataset_experiment_us_eu$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_eu)

### Control vs EU: manipulation

dataset_experiment_control_eu <- experiment_manipulation_check %>% 
  filter(experimental_groups_a=="Control" | 
           experimental_groups_a=="Blame EU")

dataset_experiment_control_eu$experimental_groups_a <- as.numeric(dataset_experiment_control_eu$experimental_groups_a)
dataset_experiment_control_eu$outcome_blame <- as.numeric(dataset_experiment_control_eu$outcome_blame)

wilcox.test(outcome_blame ~ experimental_groups_a, 
            data = dataset_experiment_control_eu) 

dataset_experiment_us_eu$experimental_groups_a <- as.numeric(dataset_experiment_us_eu$experimental_groups_a) 

cohen.d(outcome_blame ~ experimental_groups_a,
        data = dataset_experiment_us_eu)

### END SCRIPT ### END SCRIPT ### END SCRIPT ### END SCRIPT ### END SCRIPT ###



